function [stable_school_fd, stable_teach_fd, stable_school_fdva, stable_teach_fdva] = func_CF_pareto( IJipref,IJoutput, IJVAma,IJclassSize,IJVAm1, IJVAm2,IJVAmean, IJfracDisadv, IJitidx,IJjtidx,IJappyear,thetahat,yearselect, mult_fd, mult_fdva)


    %% go down to IJappyear==yearselect
    
    IJipref=IJipref(IJappyear==yearselect);
    IJitidx2015=IJitidx(IJappyear==yearselect);
    IJjtidx2015=IJjtidx(IJappyear==yearselect);
    IJVAma=IJVAma(IJappyear==yearselect);
    IJclassSize=IJclassSize(IJappyear==yearselect);
    IJoutput=IJoutput(IJappyear==yearselect);
    
    
    IJVAm1=IJVAm1(IJappyear==yearselect);
    IJVAm2=IJVAm2(IJappyear==yearselect);
    IJfracDisadv=IJfracDisadv(IJappyear==yearselect);
    IJVAmean=IJVAmean(IJappyear==yearselect);
    
    % redo indexes so that they start at 1 
    
    IJitidx=IJitidx2015-min(IJitidx2015)+1;
    IJjtidx=IJjtidx2015-min(IJjtidx2015)+1;
   
        % make balanced
         imax=max(IJitidx);       
        jlist=unique(IJjtidx);
        jkeep=randperm(max(jlist),imax)';  % select a number of classrooms that equals the number of teachers
        jselect=ismember(IJjtidx,jkeep);  % generate a logical vector telling us which to keep
        
        % now apply the selector
        
        IJipref=IJipref( jselect);
        IJVAma=IJVAma( jselect);
        IJclassSize=IJclassSize( jselect);
        IJoutput=IJoutput( jselect);
        IJVAmean=IJVAmean(jselect);
        IJVAm1=IJVAm1( jselect);
        IJVAm2=IJVAm2( jselect);
        IJfracDisadv= IJfracDisadv( jselect);
        IJitidx=IJitidx( jselect);
        IJjtidx= IJjtidx( jselect);
        
        % sort j in sequential order (to renumber), then reorder everything
        % else
        [b1,Isorter]=sort(IJjtidx); % sort by j
        inbetween=repmat((1:imax),imax,1);
        IJjtidx=inbetween(:);
        
            IJipref=IJipref(Isorter);
        IJVAma=IJVAma( Isorter);
        IJclassSize=IJclassSize( Isorter);
        IJoutput=IJoutput( Isorter);
        IJVAmean=IJVAmean(Isorter);
        IJVAm1=IJVAm1( Isorter);
        IJVAm2=IJVAm2(Isorter);
        IJfracDisadv= IJfracDisadv(Isorter);
        IJitidx=IJitidx(Isorter);
   
    
    
%% make the teacher preferences rectangular, and shocks, and then rank
    
    error_draws=evrnd(0,1,length(IJitidx),1);
    IJipref(:,1)=IJipref(:,1)+error_draws(:,1);
  

            for i=1:length(IJitidx)
                 IJipref_sq(IJitidx(i,1),IJjtidx(i,1))=IJipref(i,1); % add the type I errors
            end
       
    [B1, IJirank]=sort(IJipref_sq,2,'descend');  % this is sorted in descending order; so the first column in I is the most preferred of the teacher in the relevant row 
    
    % make the school output vector rectangular, and then rank
    % note: this breaks ties in an arbitrary way (are there ties?)
   
   for i=1:length(IJitidx)
         IJjpref_sq(IJjtidx(i,1),IJitidx(i,1))=IJoutput(i,1);    % these are J X I (useful for DA)
         IJVAma_sq(IJjtidx(i,1),IJitidx(i,1))=IJVAma(i,1); % these are J X I (useful for DA)
         IJfracDisadv_ij(IJjtidx(i,1),IJitidx(i,1))=IJfracDisadv(i,1);  % these are J x 1 (useful for first best)    
         IJijpref_sq(IJitidx(i,1),IJjtidx(i,1))=IJoutput(i,1);    % these are I X J (useful for first best)
         IJijVAma_sq(IJitidx(i,1),IJjtidx(i,1))=IJVAma(i,1);  % these are I X J (useful for first best)  
         IJclassSize_sq(IJitidx(i,1),IJjtidx(i,1))=IJclassSize(i,1);  % these are I x J (useful for first best)
         IJVAm1_sq(IJitidx(i,1),IJjtidx(i,1))=IJVAm1(i,1);  % these are I X J (useful for first best)
         IJVAm2_sq(IJitidx(i,1),IJjtidx(i,1))=IJVAm2(i,1);  % these are I X J (useful for first best)
         IJVAmean_sq(IJitidx(i,1),IJjtidx(i,1))=IJVAmean(i,1); % I x J 
         IJfracDisadv_sq(IJitidx(i,1),IJjtidx(i,1))=IJfracDisadv(i,1);  % these are I X J (useful for first best)    
      
   end
   
   IJtype1num=(1-IJfracDisadv_sq).*IJclassSize_sq;
   IJtype2num=IJfracDisadv_sq.*IJclassSize_sq;
       
    [B1, IJjrank]=sort(IJjpref_sq,2,'descend');  % this is sorted in descending order; so the first column in I is the most preferred of the school in the relevant row 
    
    % run DA
    
   [iass_jprop]= galeshapley_gender_unbalanced(IJjrank,IJirank);  % school propose
          
   [jass_iprop]= galeshapley_gender_unbalanced(IJirank,IJjrank);  % teach propose
   
   imax=max(IJitidx);
   jmax=max(IJjtidx);
    
   % figure out iass_iprop
   istuff= [ (1:length(jass_iprop))' jass_iprop];
   istuff2=sortrows(istuff,2);
      
   iass_iprop=[ zeros(imax,1)];
   for ii=1:imax
       for i=1:jmax
         if (istuff2(i,2)==ii)
            iass_iprop(ii,1)=istuff2(i,1);
         end
       end
   end     
  
  % test if the school propose and teacher propose are the same
  
  teach_school_same=min((iass_iprop==iass_jprop));

   % compute per teacher utility in teacher propose and school propose

   for i=1:min(max(IJitidx),max(IJjtidx))
     i_util_iprop(i,1)=IJipref_sq(i,iass_iprop(i,1));
   end
   
   for i=1:min(max(IJitidx),max(IJjtidx))
     i_util_jprop(i,1)=IJipref_sq(i,iass_jprop(i,1));
   end

   
   % restrict to the set of schools that appear in the stable allocation
   % (so that the bonuses don't change the set of matched schools)
   
   
  IJipref_sq=IJipref_sq(:,iass_iprop);
  IJVAm1_sq=IJVAm1_sq(:,iass_iprop);
  IJVAm2_sq=IJVAm2_sq(:,iass_iprop);
  IJtype1num=IJtype1num(:,iass_iprop); 
  IJtype2num=IJtype2num(:,iass_iprop); 
  IJclassSize_sq=IJclassSize_sq(:,iass_iprop);  
  IJfracDisadv_sq=IJfracDisadv_sq(:,iass_iprop);
  IJijVAma_sq=IJijVAma_sq(:,iass_iprop);
  IJVAmean_sq=IJVAmean_sq(:,iass_iprop);
  IJijpref_sq=IJijpref_sq(:,iass_iprop);
  
  IJjrank=IJjrank(iass_iprop,:);
    

%% Counterfactual 1: subsidize N disadvantage 

  IJipref_sq_cf=IJipref_sq+IJfracDisadv_sq.*IJclassSize_sq.* abs(thetahat(2))*mult_fd;

 [stable_teach_fd stable_school_fd] = func_CF_stable(IJipref_sq_cf,IJipref_sq_cf,IJipref_sq, IJjrank,thetahat, IJVAm1_sq, IJtype1num,  IJVAm2_sq, IJtype2num, IJijVAma_sq,IJclassSize_sq,IJitidx,IJitidx,   i_util_iprop, i_util_jprop );

%% Counterfactual 2: subsidize N disadvantage times absolute advantage

  IJipref_sq_cf=IJipref_sq+IJfracDisadv_sq.*IJclassSize_sq.*(IJVAmean_sq-min(IJVAmean)).* abs(thetahat(2))*mult_fdva;

 [stable_teach_fdva stable_school_fdva] = func_CF_stable(IJipref_sq_cf,IJipref_sq_cf,IJipref_sq, IJjrank,thetahat, IJVAm1_sq, IJtype1num,  IJVAm2_sq, IJtype2num, IJijVAma_sq,IJclassSize_sq,IJitidx,IJitidx,   i_util_iprop, i_util_jprop );

   
  end      